Scenario

load libraries

library(tidyverse)
library(readxl) 
library(grid)
library(lubridate)
library(janitor)
library(brms)
library(rstan)
library(reshape2)
library(modelr)
library(emmeans)
library(bayestestR) 
library(tidybayes) 
library(patchwork) 
library(ggpubr)
library(cowplot)
library(HDInterval)
library(ggridges) 
library(DHARMa)
library(rphylopic)
library(ggeffects)

Data import

orangemitedata_2020 <- read_delim("data/orangemitedata_2020.txt", 
                                  "\t", 
                                  escape_double = FALSE, 
                                  col_types = cols(Date = col_datetime(format = "%d/%m/%Y"),
                                                   Trapline = col_factor(levels = c("1","2", "3", "4","5", "7", "8","9")), 
                                                   Species = col_factor(levels = c("1","2", "3")), 
                                                   Sex = col_factor(levels = c("1","2")), 
                                                   Year = col_factor(levels = c("1","2", "3", "4"))), 
                                  trim_ws = TRUE)|>  
  rename(OrangeMites =`Orange Mites (Y/N)`)

smallmammaldata3 <- read_delim("data/smallmammaldata_v3.txt", 
                               "\t", 
                               escape_double = FALSE, 
                               col_types = cols(Date = col_date(format = "%Y-%m-%d"),
                                                Trapline = col_factor(levels = c("1","2", "3", "4", "5", "6", "7","8","9")),
                                                Species = col_factor(levels = c("1","2", "3")), 
                                                Sex = col_factor(levels = c("1","2")),
                                                Year = col_factor(levels = c("1","2", "3")), 
                                                X14 = col_skip(), 
                                                X15 = col_skip()), 
                               trim_ws = TRUE)|>  
  rename(OrangeMites =`Orange Mites (Y/N)`)

abund <- read_excel("data/smallmammaldata2016-2018_all.xlsx") 
abund_2020 <- read_delim("data/smallmammaldata2020_all.txt", 
                         "\t", escape_double = FALSE, col_types = cols(Date = col_date(format = "%d/%m/%Y")), 
                         trim_ws = TRUE) 

Data manipulation

orangemitedata_v2 <- bind_rows(smallmammaldata3, orangemitedata_2020) |> 
  mutate(Trapline =factor(Trapline, levels =c('4','1','2','3','5','6','7','8','9')), 
         Species =factor(Species, levels =c('2','1','3')), 
         OrangeMites =factor(OrangeMites), 
         Year =factor(Year), 
         TagLeft =case_when(TagLeft == NA ~ TagRight, 
                            TRUE ~ TagLeft), 
         Date.Julian =as.numeric(format(as.POSIXct(Date, formate="%Y-%m-%d"), "%j"))) |> 
  select(c(Date:Year, Date.Julian)) |> 
  clean_names(case="all_caps") 

abund_1 <- rbind(abund,abund_2020) |> 
  clean_names() |>
  mutate(orange_mites_y_n = toupper(orange_mites_y_n), 
         species = toupper(species), 
         orange_mites_y_n = replace_na(orange_mites_y_n, "N"), 
         year = year(date), 
         month = month(date),  
         tag_left = coalesce(tag_left,tag_right), 
         tag_left = ifelse(tag_left=="RIP", tag_right, tag_left)) |> 
  group_by(year) |> 
  mutate(first_om_appearence = min(date[orange_mites_y_n=="Y"])) |>
  ungroup() |> 
  subset(date > first_om_appearence) |> 
  mutate(month_name = month.abb[month],
         species = case_when((species == "SM") ~ "DM",
                             (species == "NFLYER") ~ "FLYER",
                             TRUE ~ species),
         trapline_group = as.numeric(trapline), #make a separate group that groups traplines based on habitat type
         trapline_group = case_when((trapline >100 & trapline<200) ~ 100,
                                    (trapline >200 & trapline<300) ~ 200,
                                    (trapline >300 & trapline<400) ~ 300,
                                    (trapline >400 & trapline<500) ~ 400,
                                    (trapline >500 & trapline<600) ~ 500,
                                    (trapline >700 & trapline<800) ~ 700,
                                    (trapline >800 & trapline<900) ~ 800,
                                    TRUE ~ 900)) |> #place new column next to the 'trapline' column to make sure it worked
  relocate(trapline_group, .after=trapline) |> 
  unite("group_ID", year, month_name, trapline_group, remove = F)

abund_2 <-  abund_1 |> 
  select(c("group_ID":"species", "tag_left", "year":"month_name")) |>
  group_by(year) |> 
  distinct(tag_left, .keep_all = T) |>
  ungroup() |>
  reshape2::dcast(group_ID~species, length) |> 
  select(-c(`TRAP MALFUNCTION`, SOREX)) |> #`TRAP DISTURBANCE`
  mutate(total =rowSums(across(-group_ID)), 
         perc_rbv = RBV/total) 
## Using month_name as value column: use value.var to override.
#join 'abund_2' dataframe to 'orangemitedata_v3' dataframe by the group_ID column 
#so that the percentage of voles within a community (per month) can be recorded for each observation
orangemitedata_v3 <- orangemitedata_v2 |> 
  mutate(year = year(DATE), 
         month = month(DATE),
         month_name = month.abb[month], 
         trapline_group= as.numeric(as.character(TRAPLINE))*100) |> 
  unite("group_ID", year,month_name,trapline_group, remove=F) |>
  left_join(abund_2[,c(1,11)], by="group_ID") |> 
  clean_names(case="all_caps") 

orangemitedata_y1y2 <- orangemitedata_v3|> 
  select(c(TRAPLINE, SPECIES, TAG_LEFT,ORANGE_MITES,YEAR,DATE_JULIAN)) |> 
  filter(TRAPLINE != 9) |>
  filter(YEAR ==1 | YEAR ==2) |>
  mutate(TRAPLINE =factor(TRAPLINE, levels=c('4','1','2','3','5','7','8')), 
         YEAR =factor(YEAR, levels =c('1','2'))) 

orangemitedata_y3y4 <- orangemitedata_v3|> 
  select(c(TRAPLINE, SPECIES, TAG_LEFT,ORANGE_MITES,YEAR,DATE_JULIAN)) |> 
  filter(TRAPLINE != 9) |> 
  filter(YEAR ==3 | YEAR ==4) |>
  mutate(TRAPLINE =factor(TRAPLINE, levels=c('4','1','2','3','5','7','8')), 
         YEAR =factor(YEAR, levels =c('3','4')))  

Descriptive statistics

temp1 <- separate(abund_2, group_ID, c("year","month","trapline.no"))|> 
  filter(year=="2016"|year=="2017")
summarise_all(temp1[4:12],funs(sum))/sum(temp1[12]) 
## Warning: `funs()` was deprecated in dplyr 0.8.0.
## ℹ Please use a list of either functions or lambdas:
## 
## # Simple named list: list(mean = mean, median = median)
## 
## # Auto named with `tibble::lst()`: tibble::lst(mean, median)
## 
## # Using lambdas list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
temp1 |> select("trapline.no","DM","RBV","WJM") |> 
  mutate(trapline.no =factor(trapline.no)) |>
  group_by(trapline.no) |> 
  summarise_each(list(sum))
## Warning: `summarise_each()` was deprecated in dplyr 0.7.0.
## ℹ Please use `across()` instead.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
temp2 <- separate(abund_2, group_ID, c("year","month","trapline.no"))|> 
  filter(year=="2018"|year=="2020")
summarise_all(temp2[4:12],funs(sum))/sum(temp2[12])
## Warning: `funs()` was deprecated in dplyr 0.8.0.
## ℹ Please use a list of either functions or lambdas:
## 
## # Simple named list: list(mean = mean, median = median)
## 
## # Auto named with `tibble::lst()`: tibble::lst(mean, median)
## 
## # Using lambdas list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
temp2 |> select("trapline.no","DM","RBV","WJM") |> 
  mutate(trapline.no =factor(trapline.no)) |>
  group_by(trapline.no) |> 
  summarise_each(list(sum))
## Warning: `summarise_each()` was deprecated in dplyr 0.7.0.
## ℹ Please use `across()` instead.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Data analysis

High-RBV year

Fit models

priors

get_priors_y1y2 <- get_prior(ORANGE_MITES~SPECIES+TRAPLINE+YEAR+scale(DATE_JULIAN, center =TRUE, scale =TRUE)+SPECIES:TRAPLINE+SPECIES:scale(DATE_JULIAN, center =TRUE, scale =TRUE)+(1|TAG_LEFT), 
                         data=orangemitedata_y1y2,
                         family=bernoulli()) ;get_priors_y1y2 
priors_y1y2 <- c(set_prior("normal(0,5)", class ="Intercept"), 
                 set_prior("normal(0,5)", class="b")) 

model formula

comm.comp1.formula <- bf(ORANGE_MITES~SPECIES+TRAPLINE+YEAR+scale(DATE_JULIAN, center=TRUE, scale =TRUE)+SPECIES:TRAPLINE+SPECIES:scale(DATE_JULIAN, center =TRUE, scale =TRUE)+(1|TAG_LEFT), 
                         family=bernoulli())

Fit model

comm.comp1.2025 <- brm(comm.comp1.formula, 
                       data =orangemitedata_y1y2, 
                       family=bernoulli(), 
                       prior = priors_y1y2,
                       iter =2000,
                       warmup = 1000, 
                       seed=123,
                       thin=2,
                       chains = 4,
                       cores = 2, 
                       save_pars = save_pars(all = TRUE), 
                       sample_prior = "yes")
saveRDS(comm.comp1.2025, "models/comm_comp1_2025.RDS")

model validations (convergence)

stan_trace(comm.comp1.2025$fit)
## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(comm.comp1.2025$fit) 
## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(comm.comp1.2025$fit) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(comm.comp1.2025$fit)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(comm.comp1.2025$fit, separate_chains = TRUE) 
## 'pars' not specified. Showing first 10 parameters by default.

model validation (DHARMa)

# not sure how...

model investigation

summary

summary(comm.comp1.2025)
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: ORANGE_MITES ~ SPECIES + TRAPLINE + YEAR + scale(DATE_JULIAN, center = TRUE, scale = TRUE) + SPECIES:TRAPLINE + SPECIES:scale(DATE_JULIAN, center = TRUE, scale = TRUE) + (1 | TAG_LEFT) 
##    Data: orangemitedata_y1y2 (Number of observations: 1446) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 2;
##          total post-warmup draws = 2000
## 
## Group-Level Effects: 
## ~TAG_LEFT (Number of levels: 603) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     2.74      0.30     2.19     3.35 1.00     1038     1600
## 
## Population-Level Effects: 
##                                                  Estimate Est.Error l-95% CI
## Intercept                                            7.04      1.04     5.16
## SPECIES1                                            -4.87      1.16    -7.22
## SPECIES3                                            -4.15      1.33    -6.75
## TRAPLINE1                                           -6.39      1.80    -9.95
## TRAPLINE2                                           -2.21      1.54    -5.13
## TRAPLINE3                                           -0.79      1.39    -3.47
## TRAPLINE5                                           -3.35      1.44    -6.24
## TRAPLINE7                                           -0.65      1.42    -3.38
## TRAPLINE8                                           -5.54      1.19    -7.95
## YEAR2                                               -0.05      0.38    -0.80
## scaleDATE_JULIANcenterEQTRUEscaleEQTRUE              1.80      0.42     1.01
## SPECIES1:TRAPLINE1                                   2.27      1.91    -1.35
## SPECIES3:TRAPLINE1                                   3.28      2.08    -0.73
## SPECIES1:TRAPLINE2                                  -0.71      1.70    -4.07
## SPECIES3:TRAPLINE2                                  -0.32      2.02    -4.28
## SPECIES1:TRAPLINE3                                  -0.70      1.53    -3.80
## SPECIES3:TRAPLINE3                                  -0.35      2.12    -4.59
## SPECIES1:TRAPLINE5                                  -0.79      1.66    -4.16
## SPECIES3:TRAPLINE5                                  -1.56      1.95    -5.47
## SPECIES1:TRAPLINE7                                  -1.62      1.57    -4.74
## SPECIES3:TRAPLINE7                                  -2.32      1.97    -6.40
## SPECIES1:TRAPLINE8                                   1.43      1.44    -1.50
## SPECIES3:TRAPLINE8                                   3.35      1.68     0.12
## SPECIES1:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE     0.05      0.45    -0.83
## SPECIES3:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE     0.27      0.54    -0.82
##                                                  u-95% CI Rhat Bulk_ESS
## Intercept                                            9.23 1.00     1201
## SPECIES1                                            -2.65 1.00     1214
## SPECIES3                                            -1.64 1.00     1450
## TRAPLINE1                                           -3.00 1.00     1774
## TRAPLINE2                                            0.84 1.00     1479
## TRAPLINE3                                            2.00 1.00     1571
## TRAPLINE5                                           -0.54 1.00     1627
## TRAPLINE7                                            2.13 1.00     1473
## TRAPLINE8                                           -3.29 1.00     1335
## YEAR2                                                0.69 1.00     1817
## scaleDATE_JULIANcenterEQTRUEscaleEQTRUE              2.66 1.00     1573
## SPECIES1:TRAPLINE1                                   6.05 1.00     1481
## SPECIES3:TRAPLINE1                                   7.30 1.00     1644
## SPECIES1:TRAPLINE2                                   2.49 1.00     1551
## SPECIES3:TRAPLINE2                                   3.58 1.00     1546
## SPECIES1:TRAPLINE3                                   2.26 1.00     1312
## SPECIES3:TRAPLINE3                                   3.70 1.00     1761
## SPECIES1:TRAPLINE5                                   2.52 1.00     1590
## SPECIES3:TRAPLINE5                                   2.33 1.00     1713
## SPECIES1:TRAPLINE7                                   1.43 1.00     1585
## SPECIES3:TRAPLINE7                                   1.40 1.00     1657
## SPECIES1:TRAPLINE8                                   4.20 1.00     1507
## SPECIES3:TRAPLINE8                                   6.62 1.00     1474
## SPECIES1:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE     0.92 1.00     1521
## SPECIES3:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE     1.34 1.00     1681
##                                                  Tail_ESS
## Intercept                                            1620
## SPECIES1                                             1426
## SPECIES3                                             1749
## TRAPLINE1                                            1684
## TRAPLINE2                                            1617
## TRAPLINE3                                            1775
## TRAPLINE5                                            1722
## TRAPLINE7                                            1787
## TRAPLINE8                                            1398
## YEAR2                                                1466
## scaleDATE_JULIANcenterEQTRUEscaleEQTRUE              1898
## SPECIES1:TRAPLINE1                                   1538
## SPECIES3:TRAPLINE1                                   1655
## SPECIES1:TRAPLINE2                                   1821
## SPECIES3:TRAPLINE2                                   1718
## SPECIES1:TRAPLINE3                                   1673
## SPECIES3:TRAPLINE3                                   1873
## SPECIES1:TRAPLINE5                                   1505
## SPECIES3:TRAPLINE5                                   1714
## SPECIES1:TRAPLINE7                                   1853
## SPECIES3:TRAPLINE7                                   1682
## SPECIES1:TRAPLINE8                                   1659
## SPECIES3:TRAPLINE8                                   1810
## SPECIES1:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE     1844
## SPECIES3:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE     1895
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

joint tests

comm.comp1.2025 |> emmeans(~ SPECIES | TRAPLINE, type ="response") |> regrid() |> joint_tests()

R-squared

comm.comp1.2025 |> bayes_R2(summary =FALSE) |> median_hdci()

post-hoc analysis

comm.comp1.2025 |> emmeans(~ SPECIES | TRAPLINE, type ="response") |> regrid() |> pairs(by="TRAPLINE") |> summary()
comm.comp1.2025 |> emmeans(~ SPECIES | TRAPLINE, type ="response") |> regrid() |> pairs(by="SPECIES") |> summary()

probabilities

comm.comp1.2025 |> emmeans(~SPECIES*TRAPLINE, type ="response")
##  SPECIES TRAPLINE response lower.HPD upper.HPD
##  2       4          0.9991  0.995720    1.0000
##  1       4          0.8955  0.682984    0.9875
##  3       4          0.9449  0.732562    0.9984
##  2       1          0.6473  0.100637    0.9995
##  1       1          0.1224  0.037927    0.2236
##  3       1          0.4401  0.133342    0.7608
##  2       2          0.9912  0.939314    1.0000
##  1       2          0.3215  0.133179    0.5222
##  3       2          0.5890  0.129636    0.9408
##  2       3          0.9978  0.987916    1.0000
##  1       3          0.6614  0.413431    0.8758
##  3       3          0.8441  0.313672    0.9999
##  2       5          0.9729  0.841554    0.9998
##  1       5          0.1232  0.024561    0.2876
##  3       5          0.1153  0.003966    0.4342
##  2       7          0.9982  0.988212    1.0000
##  1       7          0.4710  0.249584    0.7214
##  3       7          0.4781  0.090867    0.8773
##  2       8          0.8112  0.545423    0.9767
##  1       8          0.1238  0.018941    0.3104
##  3       8          0.6581  0.323746    0.9399
## 
## Results are averaged over the levels of: YEAR 
## Point estimate displayed: median 
## Results are back-transformed from the logit scale 
## HPD interval probability: 0.95

prob. comparisons

mtsqst1 <- comm.comp1.2025 |> emmeans(pairwise~SPECIES*TRAPLINE)
mtsqrt2 <- mtsqst1$contrasts |> gather_emmeans_draws() |> 
  filter(contrast %in% c(
    "SPECIES2 TRAPLINE4 - SPECIES1 TRAPLINE4", 
    "SPECIES2 TRAPLINE4 - SPECIES3 TRAPLINE4", 
    "SPECIES2 TRAPLINE2 - SPECIES1 TRAPLINE2", 
    "SPECIES2 TRAPLINE2 - SPECIES3 TRAPLINE2", 
    "SPECIES2 TRAPLINE3 - SPECIES1 TRAPLINE3", 
    "SPECIES2 TRAPLINE3 - SPECIES3 TRAPLINE3",
    "SPECIES2 TRAPLINE5 - SPECIES1 TRAPLINE5", 
    "SPECIES2 TRAPLINE5 - SPECIES3 TRAPLINE5",
    "SPECIES2 TRAPLINE7 - SPECIES1 TRAPLINE7", 
    "SPECIES2 TRAPLINE7 - SPECIES3 TRAPLINE7",
    "SPECIES2 TRAPLINE8 - SPECIES1 TRAPLINE8", 
    "SPECIES2 TRAPLINE8 - SPECIES3 TRAPLINE8",
    "SPECIES2 TRAPLINE1 - SPECIES1 TRAPLINE1", 
    "SPECIES2 TRAPLINE1 - SPECIES3 TRAPLINE1"
    ))
mtsqrt2 %>% group_by(contrast) %>% dplyr::summarise(Prob = sum(.value>0)/n()) 

partial effects plots

comm.comp1.2025 |> conditional_effects(spaghetti = TRUE, ndraws =200)

summary figure

var <- get_variables(comm.comp1.2025)
vari <- get_variables(comm.comp1.2025)[c(1,4:9,12:23)]

comm.comp1_int_draws2 <- comm.comp1.2025 |> spread_draws(!!!syms(vari)) 

comm.comp1_draws <- comm.comp1.2025 |> gather_draws(b_Intercept, b_SPECIES1, b_SPECIES3) |> 
  left_join(comm.comp1_int_draws2, by = c(".chain",".iteration",".draw")) |> 
  mutate(TRAPLINE4_mean = case_when(`.variable` == 'b_Intercept' ~ 
                                      `.value`,
                                    `.variable` != 'b_Intercept' ~ 
                                      `.value` + b_Intercept), 
         TRAPLINE1_mean = case_when(`.variable` == 'b_Intercept' ~ 
                                      `.value` + b_TRAPLINE1, 
                                    `.variable` == 'b_SPECIES1' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE1 + `b_SPECIES1:TRAPLINE1`, 
                                    `.variable` == 'b_SPECIES3' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE1 + `b_SPECIES3:TRAPLINE1`), 
         
         TRAPLINE2_mean = case_when(`.variable` == 'b_Intercept' ~ 
                                      `.value` + b_TRAPLINE2, 
                                    `.variable` == 'b_SPECIES1' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE2 + `b_SPECIES1:TRAPLINE2`, 
                                    `.variable` == 'b_SPECIES3' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE2 + `b_SPECIES3:TRAPLINE2`), 
         
         TRAPLINE3_mean = case_when(`.variable` == 'b_Intercept' ~ 
                                      `.value` + b_TRAPLINE3, 
                                    `.variable` == 'b_SPECIES1' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE3 + `b_SPECIES1:TRAPLINE3`, 
                                    `.variable` == 'b_SPECIES3' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE3 + `b_SPECIES3:TRAPLINE3`), 
         
         TRAPLINE5_mean = case_when(`.variable` == 'b_Intercept' ~ 
                                      `.value` + b_TRAPLINE5, 
                                    `.variable` == 'b_SPECIES1' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE5 + `b_SPECIES1:TRAPLINE5`, 
                                    `.variable` == 'b_SPECIES3' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE5 + `b_SPECIES3:TRAPLINE5`), 
         
         TRAPLINE7_mean = case_when(`.variable` == 'b_Intercept' ~ 
                                      `.value` + b_TRAPLINE7, 
                                    `.variable` == 'b_SPECIES1' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE7 + `b_SPECIES1:TRAPLINE7`, 
                                    `.variable` == 'b_SPECIES3' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE7 + `b_SPECIES3:TRAPLINE7`), 
         
         TRAPLINE8_mean = case_when(`.variable` == 'b_Intercept' ~ 
                                      `.value` + b_TRAPLINE8, 
                                    `.variable` == 'b_SPECIES1' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE8 + `b_SPECIES1:TRAPLINE8`, 
                                    `.variable` == 'b_SPECIES3' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE8 + `b_SPECIES3:TRAPLINE8`), 
         
         transformed_trapline4 = exp(TRAPLINE4_mean)/(1+exp(TRAPLINE4_mean)), 
         transformed_trapline1 = exp(TRAPLINE1_mean)/(1+exp(TRAPLINE1_mean)), 
         transformed_trapline2 = exp(TRAPLINE2_mean)/(1+exp(TRAPLINE2_mean)), 
         transformed_trapline3 = exp(TRAPLINE3_mean)/(1+exp(TRAPLINE3_mean)), 
         transformed_trapline5 = exp(TRAPLINE5_mean)/(1+exp(TRAPLINE5_mean)), 
         transformed_trapline7 = exp(TRAPLINE7_mean)/(1+exp(TRAPLINE7_mean)), 
         transformed_trapline8 = exp(TRAPLINE8_mean)/(1+exp(TRAPLINE8_mean))) 

comm.comp1_draws_plotting <-  comm.comp1_draws |> 
  pivot_longer(cols = starts_with("transformed"), 
               names_to = "TRAPLINE", 
               values_to = "infect_prob") %>% 
  transmute(species = case_when(`.variable` == "b_Intercept" ~ "RBV", 
                                `.variable` == "b_SPECIES1" ~ "DM", 
                                `.variable` == "b_SPECIES3" ~ "WJM"), 
            trapline = case_when(TRAPLINE == 'transformed_trapline4' ~ "trapline4", 
                                 TRAPLINE == 'transformed_trapline1' ~ "trapline1", 
                                 TRAPLINE == 'transformed_trapline2' ~ "trapline2",
                                 TRAPLINE == 'transformed_trapline3' ~ "trapline3",
                                 TRAPLINE == 'transformed_trapline5' ~ "trapline5",
                                 TRAPLINE == 'transformed_trapline7' ~ "trapline7",
                                 TRAPLINE == 'transformed_trapline8' ~ "trapline8"), 
            infect_prob = infect_prob, 
            speciesb = species, 
            traplineb = trapline, 
            chain = `.chain`, 
            iteration = `.iteration`, 
            draw_num = `.draw`) %>% 
  unite("id",speciesb,traplineb,sep = "_") %>% 
  mutate(id = factor(id, levels = c("RBV_trapline4","DM_trapline4","WJM_trapline4", 
                                    "RBV_trapline1","DM_trapline1","WJM_trapline1", 
                                    "RBV_trapline2","DM_trapline2","WJM_trapline2", 
                                    "RBV_trapline3","DM_trapline3","WJM_trapline3",
                                    "RBV_trapline5","DM_trapline5","WJM_trapline5",
                                    "RBV_trapline7","DM_trapline7","WJM_trapline7",
                                    "RBV_trapline8","DM_trapline8","WJM_trapline8"))) 

trapline_names <- c("Sugar maple hardwood", " Cut-over mixed-wood", "Dense mixed-wood",
                    "Conifer", "White pine/white spruce", "Black spruce/aspen", 
                    "White/red pine")
trapline_index <- c("trapline1","trapline2","trapline3","trapline4","trapline5",
                    "trapline7","trapline8") 
species_labels <- c("DM" = "Deer mouse","RBV" = "Red-Backed Vole","WJM" = "Woodland Jumping Mouse") 
ridge_plot1 <- comm.comp1_draws_plotting |> 
  ggplot(aes(infect_prob, trapline, fill = species)) +
  scale_fill_manual(values = alpha(c("#7C8081","#7E4640",  "#ABA159"),0.7)) +
  geom_density_ridges_gradient(scale=2, rel_min_height = 0.01) +
  facet_wrap(~species, labeller = labeller(species = species_labels))+
  theme_bw()+ 
  scale_x_continuous(breaks = seq(0, 1, by = .2))+
  xlab("Probability of infection")+ylab("Habitat type")+ 
  scale_y_discrete(limits =c("trapline1",
                             "trapline8",
                             "trapline5",
                             "trapline2", 
                             "trapline7",
                             "trapline3",
                             "trapline4"),
    labels = c("Sugar maple hardwood", 
               "White/red pine",
               "White pine/white spruce",
               "Cut-over mixed-wood",  
               "Black spruce/aspen",  
               "Dense mixed-wood",
               "Conifer"))+
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=25,face="bold"), 
        strip.text = element_text(size=18), 
        legend.position = "none", 
        panel.grid.major = element_blank(), panel.grid.minor = element_blank());ridge_plot1
## Picking joint bandwidth of 0.018
## Picking joint bandwidth of 0.0129
## Picking joint bandwidth of 0.0332

ggsave("outputs/2025_Community_composition_1.pdf", width =13, height =8, units ="in", dpi =300) 
ridge_plot1 
dev.off()
comm.comp1_int_draws2 |>
  pivot_longer(everything(), names_to = 'var', values_to = 'value') |>
  dplyr::filter(!var %in% c('.chain', '.draw', '.iteration')) |>
  ggplot(aes(x = value)) +
  geom_density() +
  facet_wrap(~ var, ncol = 1, scales = 'free_y') + 
  geom_vline(xintercept=0, color="red")

summary figure 2

date.plot <- orangemitedata_y1y2 |> 
  group_by(SPECIES, TRAPLINE) |> 
  modelr::data_grid(DATE_JULIAN=seq(from=120, 
                          to=max(orangemitedata_y1y2$DATE_JULIAN), 
                          length.out=100),
            YEAR=1) |> 
  add_fitted_draws(comm.comp1.2025,
                   n =100,
                   re_formula=NA) |>  
  #unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  ggplot(aes(x=DATE_JULIAN, color=SPECIES)) + 
  geom_line(aes(y=.value, group=paste(SPECIES, .draw)), alpha=0.05) + 
  geom_smooth(aes(y=.value, color=SPECIES), se=FALSE, size=1.5) + 
  theme_classic() + 
  facet_wrap(~TRAPLINE); date.plot
## Warning: `fitted_draws` and `add_fitted_draws` are deprecated as their names were confusing.
## - Use [add_]epred_draws() to get the expectation of the posterior predictive.
## - Use [add_]linpred_draws() to get the distribution of the linear predictor.
## - For example, you used [add_]fitted_draws(..., scale = "response"), which
##   means you most likely want [add_]epred_draws(...).
## NOTE: When updating to the new functions, note that the `model` parameter is now
##   named `object` and the `n` parameter is now named `ndraws`.
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

Low-RBV year

Fit models

priors

get_priors_y3y4 <- get_prior(ORANGE_MITES~SPECIES+TRAPLINE+YEAR+scale(DATE_JULIAN, center =TRUE, scale =TRUE)+SPECIES:TRAPLINE+SPECIES:scale(DATE_JULIAN, center =TRUE, scale =TRUE)+(1|TAG_LEFT), 
                         data=orangemitedata_y3y4,
                         family=bernoulli()) ;get_priors_y3y4 
priors_y3y4 <- c(set_prior("normal(0,5)", class ="Intercept"), 
                 set_prior("normal(0,5)", class="b")) 

model formula

comm.comp2.formula <- bf(ORANGE_MITES~SPECIES+TRAPLINE+YEAR+scale(DATE_JULIAN, center=TRUE, scale =TRUE)+SPECIES:TRAPLINE+SPECIES:scale(DATE_JULIAN, center =TRUE, scale =TRUE)+(1|TAG_LEFT), 
                         family=bernoulli())

Fit model

comm.comp2.2025 <- brm(comm.comp2.formula, 
                       data =orangemitedata_y3y4, 
                       family=bernoulli(), 
                       prior = priors_y3y4,
                       iter =8000,
                       warmup = 4000, 
                       seed=123,
                       thin=2,
                       chains = 4,
                       cores = 2, 
                       save_pars = save_pars(all = TRUE), 
                       sample_prior = "yes", 
                       control = list(adapt_delta=0.99, max_treedepth=15))
saveRDS(comm.comp2.2025, "models/comm_comp2_2025.RDS")

model validations (convergence)

stan_trace(comm.comp2.2025$fit)
## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(comm.comp2.2025$fit) 
## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(comm.comp2.2025$fit) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(comm.comp2.2025$fit)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(comm.comp2.2025$fit, separate_chains = TRUE) 
## 'pars' not specified. Showing first 10 parameters by default.

model investigation

summary

summary(comm.comp2.2025)
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: ORANGE_MITES ~ SPECIES + TRAPLINE + YEAR + scale(DATE_JULIAN, center = TRUE, scale = TRUE) + SPECIES:TRAPLINE + SPECIES:scale(DATE_JULIAN, center = TRUE, scale = TRUE) + (1 | TAG_LEFT) 
##    Data: orangemitedata_y3y4 (Number of observations: 1125) 
##   Draws: 4 chains, each with iter = 8000; warmup = 4000; thin = 2;
##          total post-warmup draws = 8000
## 
## Group-Level Effects: 
## ~TAG_LEFT (Number of levels: 436) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     4.46      0.69     3.27     5.96 1.00     4538     6076
## 
## Population-Level Effects: 
##                                                  Estimate Est.Error l-95% CI
## Intercept                                            4.97      1.73     1.65
## SPECIES1                                            -8.29      1.89   -12.13
## SPECIES3                                           -11.03      3.22   -17.68
## TRAPLINE1                                           -3.43      3.69   -10.71
## TRAPLINE2                                           -0.34      2.40    -5.01
## TRAPLINE3                                           -1.31      2.70    -6.58
## TRAPLINE5                                           -1.11      2.13    -5.19
## TRAPLINE7                                           -0.12      2.07    -4.25
## TRAPLINE8                                           -3.95      3.06   -10.02
## YEAR4                                               -2.84      0.84    -4.63
## scaleDATE_JULIANcenterEQTRUEscaleEQTRUE              1.08      0.59    -0.03
## SPECIES1:TRAPLINE1                                  -3.40      3.72   -10.61
## SPECIES3:TRAPLINE1                                   0.00      4.96    -9.85
## SPECIES1:TRAPLINE2                                  -0.26      2.51    -5.15
## SPECIES3:TRAPLINE2                                  -0.09      5.04    -9.87
## SPECIES1:TRAPLINE3                                  -1.58      2.84    -7.19
## SPECIES3:TRAPLINE3                                   0.04      4.99    -9.65
## SPECIES1:TRAPLINE5                                  -1.26      2.40    -6.03
## SPECIES3:TRAPLINE5                                  -2.23      4.29   -10.98
## SPECIES1:TRAPLINE7                                  -1.98      2.27    -6.61
## SPECIES3:TRAPLINE7                                  -4.47      3.77   -11.96
## SPECIES1:TRAPLINE8                                   1.10      3.18    -5.05
## SPECIES3:TRAPLINE8                                  -2.27      4.21   -10.77
## SPECIES1:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE    -0.50      0.63    -1.75
## SPECIES3:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE    -1.21      2.47    -6.11
##                                                  u-95% CI Rhat Bulk_ESS
## Intercept                                            8.48 1.00     6528
## SPECIES1                                            -4.72 1.00     6218
## SPECIES3                                            -5.03 1.00     6833
## TRAPLINE1                                            3.77 1.00     7464
## TRAPLINE2                                            4.27 1.00     6261
## TRAPLINE3                                            4.00 1.00     7257
## TRAPLINE5                                            3.02 1.00     5663
## TRAPLINE7                                            4.00 1.00     6312
## TRAPLINE8                                            1.79 1.00     7430
## YEAR4                                               -1.34 1.00     6184
## scaleDATE_JULIANcenterEQTRUEscaleEQTRUE              2.25 1.00     7013
## SPECIES1:TRAPLINE1                                   3.84 1.00     7646
## SPECIES3:TRAPLINE1                                   9.64 1.00     6523
## SPECIES1:TRAPLINE2                                   4.74 1.00     6063
## SPECIES3:TRAPLINE2                                   9.82 1.00     6092
## SPECIES1:TRAPLINE3                                   3.96 1.00     6763
## SPECIES3:TRAPLINE3                                   9.44 1.00     6022
## SPECIES1:TRAPLINE5                                   3.35 1.00     5956
## SPECIES3:TRAPLINE5                                   5.90 1.00     6866
## SPECIES1:TRAPLINE7                                   2.42 1.00     6431
## SPECIES3:TRAPLINE7                                   2.69 1.00     7081
## SPECIES1:TRAPLINE8                                   7.32 1.00     7088
## SPECIES3:TRAPLINE8                                   5.70 1.00     6517
## SPECIES1:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE     0.69 1.00     6881
## SPECIES3:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE     3.70 1.00     7080
##                                                  Tail_ESS
## Intercept                                            7175
## SPECIES1                                             6756
## SPECIES3                                             7320
## TRAPLINE1                                            6592
## TRAPLINE2                                            6549
## TRAPLINE3                                            7339
## TRAPLINE5                                            6458
## TRAPLINE7                                            7181
## TRAPLINE8                                            6645
## YEAR4                                                6984
## scaleDATE_JULIANcenterEQTRUEscaleEQTRUE              7252
## SPECIES1:TRAPLINE1                                   6686
## SPECIES3:TRAPLINE1                                   6019
## SPECIES1:TRAPLINE2                                   6949
## SPECIES3:TRAPLINE2                                   5818
## SPECIES1:TRAPLINE3                                   6885
## SPECIES3:TRAPLINE3                                   6171
## SPECIES1:TRAPLINE5                                   6406
## SPECIES3:TRAPLINE5                                   6962
## SPECIES1:TRAPLINE7                                   6738
## SPECIES3:TRAPLINE7                                   7366
## SPECIES1:TRAPLINE8                                   6872
## SPECIES3:TRAPLINE8                                   6773
## SPECIES1:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE     7250
## SPECIES3:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE     7422
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

joint tests

comm.comp2.2025 |> emmeans(~ SPECIES | TRAPLINE, type ="response") |> regrid() |> joint_tests()

R-squared

comm.comp2.2025 |> bayes_R2(summary =FALSE) |> median_hdci()

post-hoc analysis

comm.comp2.2025 |> emmeans(~ SPECIES | TRAPLINE, type ="response") |> regrid() |> pairs(by="TRAPLINE") |> summary()
comm.comp2.2025 |> emmeans(~ SPECIES | TRAPLINE, type ="response") |> regrid() |> pairs(by="SPECIES") |> summary()

probabilities

comm.comp2.2025 |> emmeans(~SPECIES*TRAPLINE, type ="response")
##  SPECIES TRAPLINE response lower.HPD upper.HPD
##  2       4        9.71e-01  7.03e-01  0.999963
##  1       4        9.05e-03  1.71e-05  0.056271
##  3       4        6.56e-04  0.00e+00  0.099588
##  2       1        5.31e-01  1.50e-03  0.999997
##  1       1        1.14e-05  0.00e+00  0.000206
##  3       1        1.83e-05  0.00e+00  0.660663
##  2       2        9.60e-01  4.14e-01  0.999997
##  1       2        5.19e-03  2.36e-05  0.023954
##  3       2        3.82e-04  0.00e+00  0.920104
##  2       3        8.99e-01  1.14e-01  0.999998
##  1       3        5.53e-04  4.00e-07  0.004829
##  3       3        1.61e-04  0.00e+00  0.858584
##  2       5        9.17e-01  4.00e-01  0.999993
##  1       5        9.24e-04  3.00e-07  0.007364
##  3       5        2.71e-05  0.00e+00  0.031493
##  2       7        9.66e-01  6.97e-01  0.999951
##  1       7        1.16e-03  2.00e-06  0.007082
##  3       7        7.60e-06  0.00e+00  0.001812
##  2       8        4.15e-01  4.70e-06  0.992098
##  1       8        5.67e-04  0.00e+00  0.005284
##  3       8        1.80e-06  0.00e+00  0.001603
## 
## Results are averaged over the levels of: YEAR 
## Point estimate displayed: median 
## Results are back-transformed from the logit scale 
## HPD interval probability: 0.95
mtsqst1 <- comm.comp2.2025 |> emmeans(pairwise~SPECIES*TRAPLINE)
mtsqrt2 <- mtsqst1$contrasts |> gather_emmeans_draws() |> 
  filter(contrast %in% c(
    "SPECIES2 TRAPLINE4 - SPECIES1 TRAPLINE4", 
    #"SPECIES2 TRAPLINE4 - SPECIES3 TRAPLINE4", 
    "SPECIES2 TRAPLINE2 - SPECIES1 TRAPLINE2", 
    #"SPECIES2 TRAPLINE2 - SPECIES3 TRAPLINE2", 
    "SPECIES2 TRAPLINE3 - SPECIES1 TRAPLINE3", 
    #"SPECIES2 TRAPLINE3 - SPECIES3 TRAPLINE3",
    #"SPECIES2 TRAPLINE5 - SPECIES1 TRAPLINE5", 
    #"SPECIES2 TRAPLINE5 - SPECIES3 TRAPLINE5",
    "SPECIES2 TRAPLINE7 - SPECIES1 TRAPLINE7", 
    #"SPECIES2 TRAPLINE7 - SPECIES3 TRAPLINE7",
    "SPECIES2 TRAPLINE8 - SPECIES1 TRAPLINE8" 
    #"SPECIES2 TRAPLINE8 - SPECIES3 TRAPLINE8",
    #"SPECIES2 TRAPLINE1 - SPECIES1 TRAPLINE1", 
    #"SPECIES2 TRAPLINE1 - SPECIES3 TRAPLINE1"
    ))
mtsqrt2 %>% group_by(contrast) %>% dplyr::summarise(Prob = sum(.value>0)/n())  

partial effects plots

comm.comp2.2025 |> conditional_effects(spaghetti = TRUE, ndraws =200)

summary figure

var <- get_variables(comm.comp2.2025)
vari <- get_variables(comm.comp2.2025)[c(1,4:9,12:23)]

comm.comp2_int_draws2 <- comm.comp2.2025 |> spread_draws(!!!syms(vari)) 

comm.comp2_draws <- comm.comp2.2025 |> gather_draws(b_Intercept, b_SPECIES1, b_SPECIES3) |> 
  left_join(comm.comp2_int_draws2, by = c(".chain",".iteration",".draw")) |> 
  mutate(TRAPLINE4_mean = case_when(`.variable` == 'b_Intercept' ~ 
                                      `.value`,
                                    `.variable` != 'b_Intercept' ~ 
                                      `.value` + b_Intercept), 
         TRAPLINE1_mean = case_when(`.variable` == 'b_Intercept' ~ 
                                      `.value` + b_TRAPLINE1, 
                                    `.variable` == 'b_SPECIES1' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE1 + `b_SPECIES1:TRAPLINE1`, 
                                    `.variable` == 'b_SPECIES3' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE1 + `b_SPECIES3:TRAPLINE1`), 
         
         TRAPLINE2_mean = case_when(`.variable` == 'b_Intercept' ~ 
                                      `.value` + b_TRAPLINE2, 
                                    `.variable` == 'b_SPECIES1' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE2 + `b_SPECIES1:TRAPLINE2`, 
                                    `.variable` == 'b_SPECIES3' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE2 + `b_SPECIES3:TRAPLINE2`), 
         
         TRAPLINE3_mean = case_when(`.variable` == 'b_Intercept' ~ 
                                      `.value` + b_TRAPLINE3, 
                                    `.variable` == 'b_SPECIES1' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE3 + `b_SPECIES1:TRAPLINE3`, 
                                    `.variable` == 'b_SPECIES3' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE3 + `b_SPECIES3:TRAPLINE3`), 
         
         TRAPLINE5_mean = case_when(`.variable` == 'b_Intercept' ~ 
                                      `.value` + b_TRAPLINE5, 
                                    `.variable` == 'b_SPECIES1' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE5 + `b_SPECIES1:TRAPLINE5`, 
                                    `.variable` == 'b_SPECIES3' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE5 + `b_SPECIES3:TRAPLINE5`), 
         
         TRAPLINE7_mean = case_when(`.variable` == 'b_Intercept' ~ 
                                      `.value` + b_TRAPLINE7, 
                                    `.variable` == 'b_SPECIES1' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE7 + `b_SPECIES1:TRAPLINE7`, 
                                    `.variable` == 'b_SPECIES3' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE7 + `b_SPECIES3:TRAPLINE7`), 
         
         TRAPLINE8_mean = case_when(`.variable` == 'b_Intercept' ~ 
                                      `.value` + b_TRAPLINE8, 
                                    `.variable` == 'b_SPECIES1' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE8 + `b_SPECIES1:TRAPLINE8`, 
                                    `.variable` == 'b_SPECIES3' ~ 
                                      `.value` + b_Intercept + 
                                      b_TRAPLINE8 + `b_SPECIES3:TRAPLINE8`), 
         
         transformed_trapline4 = exp(TRAPLINE4_mean)/(1+exp(TRAPLINE4_mean)), 
         transformed_trapline1 = exp(TRAPLINE1_mean)/(1+exp(TRAPLINE1_mean)), 
         transformed_trapline2 = exp(TRAPLINE2_mean)/(1+exp(TRAPLINE2_mean)), 
         transformed_trapline3 = exp(TRAPLINE3_mean)/(1+exp(TRAPLINE3_mean)), 
         transformed_trapline5 = exp(TRAPLINE5_mean)/(1+exp(TRAPLINE5_mean)), 
         transformed_trapline7 = exp(TRAPLINE7_mean)/(1+exp(TRAPLINE7_mean)), 
         transformed_trapline8 = exp(TRAPLINE8_mean)/(1+exp(TRAPLINE8_mean))) 

comm.comp2_draws_plotting <-  comm.comp2_draws %>% 
  pivot_longer(cols = starts_with("transformed"), 
               names_to = "TRAPLINE", 
               values_to = "infect_prob") %>% 
  transmute(species = case_when(`.variable` == "b_Intercept" ~ "RBV", 
                                `.variable` == "b_SPECIES1" ~ "DM", 
                                `.variable` == "b_SPECIES3" ~ "WJM"), 
            trapline = case_when(TRAPLINE == 'transformed_trapline4' ~ "trapline4", 
                                 TRAPLINE == 'transformed_trapline1' ~ "trapline1", 
                                 TRAPLINE == 'transformed_trapline2' ~ "trapline2",
                                 TRAPLINE == 'transformed_trapline3' ~ "trapline3",
                                 TRAPLINE == 'transformed_trapline5' ~ "trapline5",
                                 TRAPLINE == 'transformed_trapline7' ~ "trapline7",
                                 TRAPLINE == 'transformed_trapline8' ~ "trapline8"), 
            infect_prob = infect_prob, 
            speciesb = species, 
            traplineb = trapline, 
            chain = `.chain`, 
            iteration = `.iteration`, 
            draw_num = `.draw`) %>% 
  unite("id",speciesb,traplineb,sep = "_") %>% 
  mutate(id = factor(id, levels = c("RBV_trapline4","DM_trapline4","WJM_trapline4", 
                                    "RBV_trapline1","DM_trapline1","WJM_trapline1", 
                                    "RBV_trapline2","DM_trapline2","WJM_trapline2", 
                                    "RBV_trapline3","DM_trapline3","WJM_trapline3",
                                    "RBV_trapline5","DM_trapline5","WJM_trapline5",
                                    "RBV_trapline7","DM_trapline7","WJM_trapline7",
                                    "RBV_trapline8","DM_trapline8","WJM_trapline8"))) 

trapline_names <- c("Sugar maple hardwood", " Cut-over mixed-wood", "Dense mixed-wood",
                    "Conifer", "White pine/white spruce", "Black spruce/aspen", 
                    "White/red pine")
trapline_index <- c("trapline1","trapline2","trapline3","trapline4","trapline5",
                    "trapline7","trapline8") 
species_labels <- c("DM" = "Deer mouse","RBV" = "Red-Backed Vole","WJM" = "Woodland Jumping Mouse") 
ridge_plot2 <- comm.comp2_draws_plotting %>% 
  ggplot(aes(infect_prob, trapline, fill = species)) +
  scale_fill_manual(values = alpha(c("#7C8081","#7E4640",  "#ABA159"),0.7)) +
  geom_density_ridges_gradient(scale = 2, rel_min_height = 0.01)+
  facet_wrap(~species, labeller = labeller(species = species_labels))+
  theme_bw()+ 
  scale_x_continuous(breaks = seq(0, 1, by = .2))+
  xlab("P(Infection)")+ylab("Habitat type")+ 
  scale_y_discrete(limits =c("trapline1",
                             "trapline8",
                             "trapline5",
                             "trapline2", 
                             "trapline7",
                             "trapline3",
                             "trapline4"),
    labels = c("Sugar maple hardwood", 
               "White/red pine",
               "White pine/white spruce",
               "Cut-over mixed-wood",  
               "Black spruce/aspen",  
               "Dense mixed-wood",
               "Conifer"))+
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=25,face="bold"), 
        strip.text = element_text(size=18), 
        legend.position = "none", 
        panel.grid.major = element_blank(), panel.grid.minor = element_blank());ridge_plot2
## Picking joint bandwidth of 0.00185
## Picking joint bandwidth of 0.0204
## Picking joint bandwidth of 0.00299

              #img = img1, x=0.2, y=7.8, height =2, fill ="#7C8081", 
              # color = "black", 
               #alpha = 0.6)
ggsave("outputs/2025_Community_composition_2.pdf", width =13, height =8, units ="in", dpi =300) 
ridge_plot2 
dev.off()

summary figure 2

date.plot <- orangemitedata_y3y4 |> 
  group_by(SPECIES, TRAPLINE) |> 
  modelr::data_grid(DATE_JULIAN=seq(from=120, 
                          to=max(orangemitedata_y3y4$DATE_JULIAN), 
                          length.out=100),
            YEAR=3) |> 
  add_fitted_draws(comm.comp2.2025,
                   n =100,
                   re_formula=NA) |>  
  #unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  ggplot(aes(x=DATE_JULIAN, color=SPECIES)) + 
  geom_line(aes(y=.value, group=paste(SPECIES, .draw)), alpha=0.05) + 
  geom_smooth(aes(y=.value, color=SPECIES), se=FALSE, size=1.5) + 
  theme_classic() + 
  facet_wrap(~TRAPLINE); date.plot
## Warning: `fitted_draws` and `add_fitted_draws` are deprecated as their names were confusing.
## - Use [add_]epred_draws() to get the expectation of the posterior predictive.
## - Use [add_]linpred_draws() to get the distribution of the linear predictor.
## - For example, you used [add_]fitted_draws(..., scale = "response"), which
##   means you most likely want [add_]epred_draws(...).
## NOTE: When updating to the new functions, note that the `model` parameter is now
##   named `object` and the `n` parameter is now named `ndraws`.
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

RBV community presence

Fit models

Data manipulation

orangemitedata_perc.RBV <- orangemitedata_v3|> 
  select(c(TRAPLINE, SPECIES, TAG_LEFT,ORANGE_MITES,YEAR,DATE_JULIAN,PERC_RBV)) |> 
  filter(TRAPLINE != 9) |>
  mutate(TRAPLINE =factor(TRAPLINE, levels=c('4','1','2','3','5','6','7','8'))) 

priors

get_rbv.perc <- get_prior(ORANGE_MITES~PERC_RBV*SPECIES+(1|TAG_LEFT), 
                family=bernoulli(), 
                data=orangemitedata_perc.RBV);get_rbv.perc 
rbv.perc_priors <- c(set_prior("normal(0,5)", class ="Intercept"), 
                 set_prior("normal(0,5)", class="b")) 

model formula

rbv.perc_formula <- bf(ORANGE_MITES~PERC_RBV*SPECIES+(1|TAG_LEFT), 
                family=bernoulli())

Fit model

rbv.perc.2025 <- brm(rbv.perc_formula, 
                       data =orangemitedata_perc.RBV, 
                       family=bernoulli(), 
                       prior =rbv.perc_priors,
                       iter =2000,
                       warmup =1000, 
                       seed=123,
                       thin=2,
                       chains = 4,
                       cores = 2, 
                       save_pars = save_pars(all = TRUE), 
                       sample_prior = "yes", 
                       control = list(adapt_delta=0.99, max_treedepth=15))
saveRDS(rbv.perc.2025, "models/rbv_perc_2025.RDS")

model validations (convergence)

stan_trace(rbv.perc.2025$fit)
## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(rbv.perc.2025$fit) 
## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(rbv.perc.2025$fit) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(rbv.perc.2025$fit)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(rbv.perc.2025$fit, separate_chains = TRUE) 
## 'pars' not specified. Showing first 10 parameters by default.

model investigation

summary

summary(rbv.perc.2025)
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: ORANGE_MITES ~ PERC_RBV * SPECIES + (1 | TAG_LEFT) 
##    Data: orangemitedata_perc.RBV (Number of observations: 2571) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 2;
##          total post-warmup draws = 2000
## 
## Group-Level Effects: 
## ~TAG_LEFT (Number of levels: 1037) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     3.01      0.24     2.57     3.53 1.01      907     1373
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept             1.51      0.61     0.38     2.70 1.00     1190     1629
## PERC_RBV              8.24      1.72     4.98    11.66 1.00     1411     1606
## SPECIES1             -5.02      0.69    -6.44    -3.71 1.01     1068     1407
## SPECIES3             -3.85      0.78    -5.35    -2.38 1.01     1382     1828
## PERC_RBV:SPECIES1    -2.81      1.95    -6.64     0.84 1.00     1401     1441
## PERC_RBV:SPECIES3     1.33      2.33    -3.24     6.14 1.00     1478     1753
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

hpd.summary

rbv.perc.2025 |> emmeans(~SPECIES*PERC_RBV, type = "response") |> 
  hpd.summary(point.est =median)

R-squared

rbv.perc.2025 |> bayes_R2(summary =FALSE) |> median_hdci()

post-hoc analysis

rbv.perc.2025 |> emtrends(var = "PERC_RBV", type = "response") |> regrid() |> summary(by = NULL, adjust = "tukey", infer=TRUE) 
rbv.perc.2025 |> emtrends(var = "PERC_RBV", type = "response") |> 
  pairs(by = "PERC_RBV") |>
  summary(by = NULL, adjust = "tukey", infer=TRUE) 
rbv.perc.2025 |> emmeans(~SPECIES*PERC_RBV, type = "response")
##  SPECIES PERC_RBV response lower.HPD upper.HPD
##  2          0.156   0.9412    0.8866    0.9810
##  1          0.156   0.0652    0.0402    0.0961
##  3          0.156   0.3027    0.1647    0.4472
## 
## Point estimate displayed: median 
## Results are back-transformed from the logit scale 
## HPD interval probability: 0.95
rbv.perc.2025 |> emmeans(~SPECIES*PERC_RBV, type = "response") |> 
  pairs(by = "PERC_RBV") |>
  summary(by = NULL, 
          adjust = "tukey", 
          infer=TRUE)  

partial effects plots

rbv.perc.2025 |> conditional_effects(spaghetti = TRUE, ndraws =200)

Summary

perc_rbv.figure <- orangemitedata_v3 |>
  group_by(SPECIES) |> 
  modelr::data_grid(PERC_RBV=seq(from = min(orangemitedata_v3$PERC_RBV), 
                         to=max(orangemitedata_v3$PERC_RBV), 
                         length.out = 101)) |> 
  add_linpred_draws(rbv.perc.2025, ndraws=100,
                   re_formula = NA) |>
  mutate(.linpred =exp(.linpred)/(1+exp(.linpred))) |>
  ggplot(aes(x=PERC_RBV, color=SPECIES))+
  geom_line(aes(y=.linpred, group = paste(SPECIES, .draw)), alpha=0.2)+ 
  geom_smooth(aes(y=.linpred, color=SPECIES), se=F, size=1.5)+theme_classic()+ 
  xlab("% red backed vole in community composition")+ylab("Probability of infection")+ 
  scale_color_manual(labels=c("Red backed vole", "Deer mouse", "Woodland jumping mouse"), 
                     values=c("#7E4640","#7C8081","#ABA159"))   +
  theme(legend.position = "bottom", 
        legend.title = element_blank(), 
        axis.title.x = element_text(size=25), 
        axis.text.x = element_text(size=18), 
        axis.title.y = element_text(size=25), 
        axis.text.y = element_text(size=18), 
        legend.text = element_text(size=15), 
        axis.line = element_line(size=1.2)) 
## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
  #geom_vline(xintercept = 0.2, color="black", alpha=0.4, size=1); prp 
perc_rbv.figure 
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

ggsave("outputs/2025_perc_rbv_pop_model.pdf", width =13, height =8, units ="in", dpi =300) 
perc_rbv.figure 
dev.off()

Body condition

Load additional packages

library(vegan)

Data manipulation

sampling<- read_csv("data/sampling_periods.csv")
## Rows: 51 Columns: 2
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (2): julian_date, sampling_period
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
bc <- read_excel("data/Small Mammal Project 2016 (ALL)_cleaned.xlsx", 
                 sheet = "Falls Lines") |> 
  clean_names() |> #use the janitor package to clean up the names so everything's consistent
  filter(age == 'A') |> #subset to only adults JV
  filter(!(repro_cond == 'PERF/PREG' | repro_cond == 'Preg' | repro_cond == 'PREG' |
             repro_cond == 'Preg/Lac'| repro_cond == 'PREG/LAC' | 
             repro_cond == 'PREG/PERF')) |> #subset to remove preg JV
  mutate(skull_length = as.numeric(skull_length), tag_left = as.factor(tag_left)) |> #convert skull length to numeric like the rest of the variables JV 
  filter(species == 'DM' | species == 'RBV') |> #subset the df to only the species we're interested in
  mutate(julian_date = yday(date)) |> #group into sampling periods based on date
  #connect the sampling csv with this df by matching the julian dates
  #from each together
  left_join(sampling, by = "julian_date") |> 
  #so now we have a whole column of sampling periods based on the julian date
  #for each observation
  #now we order by individual
  arrange(tag_left) |> 
  #the group by individual and sampling period
  group_by(tag_left, sampling_period) |> 
  ##create average for weight JV
  mutate(avg_weight_sampling_period = mean(weight, na.rm = TRUE)) |>
  mutate(tag_left_num = as.numeric(tag_left)) |>
  distinct(avg_weight_sampling_period, tag_left_num, .keep_all = TRUE) |>
  #now remove that grouping variable so we can group by individual over the
  #whole summer now
  ungroup() |> 
  #by individual across all samples
  group_by(tag_left) |> 
  #make a new column to average measures we care about JV
  mutate(avg_tail_summer = mean(tail, na.rm = TRUE),
         avg_hindfoot_summer = mean(hindfoot, na.rm = TRUE),
         avg_ear_summer = mean(ear, na.rm = TRUE),
         avg_skull_width_summer = mean(skull_width, na.rm = TRUE),
         avg_skull_length_summer = mean(skull_length, na.rm = TRUE)) |>
  ##creating infection stages col JV
  mutate(group = case_when(any(orange_mites_y_n == "Y") ~ "Parasitized", TRUE ~ "Uninfected")) |>
  ungroup() |>
  mutate(infection_stage = case_when((group == "Parasitized" & orange_mites_y_n == "Y") ~ "After",
                                     (group == "Parasitized" & orange_mites_y_n == "N") ~ "Before",
                                     (group == "Uninfected") ~ "Before",
                                     TRUE ~ "NA")) |>
  ##There were two individuals that switched from having mites to not having mites
  ##The code below makes sure they are still listed as "after" when they lose their mites JV
  mutate(infection_stage = case_when((tag_left == 50468 & julian_date > 196) ~ "After",
                                     (tag_left == 50968 & julian_date > 207) ~ "After",
                                     TRUE ~ infection_stage)) |>
  mutate(group = as.factor(group), infection_stage = as.factor(infection_stage)) |>
  filter(infection_stage != "NA") |>
  unite(infection, c("group", "infection_stage"), remove = FALSE) |>
  mutate(infection = as.factor(infection)) |>
  ##removing NA's JV
  drop_na(avg_skull_width_summer, avg_skull_length_summer, avg_hindfoot_summer,
          avg_weight_sampling_period) |>
  ##formatting predictors JV
  mutate(sex = as.factor(sex)) 

Split dataframes

rbv = bc |> 
  filter(species == "RBV") |> 
  mutate(sex = as.factor(sex)) 

dm = bc |> 
  filter(species == "DM") |> 
  mutate(sex = as.factor(sex)) 

Red-backed voles

Princple component analysis (PCA)

pca_rbv <- princomp(~ avg_hindfoot_summer + avg_skull_length_summer + avg_skull_width_summer, 
                    rbv, 
                    cor = TRUE)
loadings(pca_rbv)
## 
## Loadings:
##                         Comp.1 Comp.2 Comp.3
## avg_hindfoot_summer      0.562  0.747  0.355
## avg_skull_length_summer  0.599        -0.798
## avg_skull_width_summer   0.570 -0.661  0.487
## 
##                Comp.1 Comp.2 Comp.3
## SS loadings     1.000  1.000  1.000
## Proportion Var  0.333  0.333  0.333
## Cumulative Var  0.333  0.667  1.000
summary(pca_rbv)
## Importance of components:
##                           Comp.1    Comp.2    Comp.3
## Standard deviation     1.3071989 0.8296390 0.7764858
## Proportion of Variance 0.5695896 0.2294336 0.2009768
## Cumulative Proportion  0.5695896 0.7990232 1.0000000
screeplot(pca_rbv)

eigenvals(pca_rbv) 
##    Comp.1    Comp.2    Comp.3 
## 1.7087689 0.6883008 0.6029303
body_size_rbv <- (pca_rbv$scores[,1])

Data exploration

plot(avg_weight_sampling_period ~ body_size_rbv, 
     data =rbv)

Extracting residuals

bodycond_correlation_rbv <- lm(avg_weight_sampling_period ~ body_size_rbv, 
                               data =rbv)
rbv$bodycondition <- bodycond_correlation_rbv$residuals
hist(rbv$bodycondition)

priors

get_RBVbodycondition_priors  <- get_prior(bodycondition ~ infection+sex*scale(julian_date)+ (1|tag_left), 
                      data=rbv,
                      family=gaussian()); get_RBVbodycondition_priors 
RBVbodycondition_priors <- c(set_prior("normal(0,5)", class ="Intercept"),
                             set_prior("normal(0,5)", class="b")) 

model formula

RBVbodycondition_formula <- bf(bodycondition ~ infection+sex*
                                 scale(julian_date, center=TRUE, scale=TRUE)+(1|tag_left), 
                               family =gaussian())

Fit model (work in progress)

RBVbodycondition <- brm(RBVbodycondition_formula, 
                       data =rbv, 
                       family=bernoulli(), 
                       prior =RBVbodycondition_priors,
                       iter =2000,
                       warmup =1000, 
                       seed=123,
                       thin=2,
                       chains = 4,
                       cores = 2, 
                       save_pars = save_pars(all = TRUE), 
                       sample_prior = "yes", 
                       control = list(adapt_delta=0.99, max_treedepth=15))
vsaveRDS(RBVbodycondition, "models/RBVbodycondition_2025.RDS")

model validations (convergence)

stan_trace(RBVbodycondition$fit)
## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(RBVbodycondition$fit) 
## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(RBVbodycondition$fit) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(RBVbodycondition$fit)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(RBVbodycondition$fit, separate_chains = TRUE) 
## 'pars' not specified. Showing first 10 parameters by default.

model investigation

summary
summary(RBVbodycondition)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: bodycondition ~ infection + sex * scale(julian_date, center = TRUE, scale = TRUE) + (1 | tag_left) 
##    Data: rbv (Number of observations: 115) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 2;
##          total post-warmup draws = 2000
## 
## Group-Level Effects: 
## ~tag_left (Number of levels: 86) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     3.18      0.37     2.49     3.93 1.01      588     1022
## 
## Population-Level Effects: 
##                                              Estimate Est.Error l-95% CI
## Intercept                                        1.08      0.68    -0.21
## infectionParasitized_Before                     -0.96      1.13    -3.14
## infectionUninfected_Before                       0.44      1.18    -1.88
## sexM                                            -2.05      0.83    -3.63
## scalejulian_datecenterEQTRUEscaleEQTRUE          0.98      0.60    -0.25
## sexM:scalejulian_datecenterEQTRUEscaleEQTRUE    -0.95      0.69    -2.31
##                                              u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept                                        2.45 1.00     1074     1506
## infectionParasitized_Before                      1.25 1.00     1581     1743
## infectionUninfected_Before                       2.68 1.00     1415     1650
## sexM                                            -0.44 1.00     1195     1465
## scalejulian_datecenterEQTRUEscaleEQTRUE          2.10 1.00     1434     1617
## sexM:scalejulian_datecenterEQTRUEscaleEQTRUE     0.37 1.00     1584     1606
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.88      0.28     1.42     2.49 1.01      553      852
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
posteriors <- as.data.frame(RBVbodycondition) 
covariate_samples <- posteriors |> select(b_Intercept, b_infectionParasitized_Before,b_infectionUninfected_Before) 
combined_effect <- rowSums(covariate_samples) 
hdi(combined_effect, credMass =0.95)
##     lower     upper 
## -2.958808  4.128201 
## attr(,"credMass")
## [1] 0.95
R-squared
RBVbodycondition |> bayes_R2(summary =FALSE) |> median_hdci()

post-hoc analysis

Infection
RBVbodycondition |> emmeans(~infection, type = "response") 
##  infection           emmean lower.HPD upper.HPD
##  Parasitized_After   0.0521    -0.827     0.962
##  Parasitized_Before -0.8898    -3.117     1.222
##  Uninfected_Before   0.5414    -1.525     2.531
## 
## Results are averaged over the levels of: sex 
## Point estimate displayed: median 
## HPD interval probability: 0.95
RBVbodycondition |> emmeans(~infection, type = "response") |> 
  pairs() |> summary(by = NULL, adjust = "tukey", infer=TRUE) 
sex
RBVbodycondition |> emmeans(~sex, type = "response")
## NOTE: Results may be misleading due to involvement in interactions
##  sex emmean lower.HPD upper.HPD
##  F    0.921    -0.708     2.482
##  M   -1.144    -2.311     0.147
## 
## Results are averaged over the levels of: infection 
## Point estimate displayed: median 
## HPD interval probability: 0.95
RBVbodycondition |> emmeans(~sex, type = "response") |> 
  pairs() |>
  summary(by = NULL, 
          adjust = "tukey", 
          infer=TRUE)  
## NOTE: Results may be misleading due to involvement in interactions
julian date
RBVbodycondition |> emtrends(var = "julian_date", type = "response") |> 
  summary(by = NULL, adjust = "tukey", infer=TRUE) 

partial effects plots

RBVbodycondition |> conditional_effects(spaghetti = TRUE, ndraws =200)

summary figure

Deer mice

Princple component analysis (PCA)

pca_dm <- princomp(~ avg_hindfoot_summer + avg_skull_length_summer + avg_skull_width_summer,
                   dm, 
                   cor = TRUE)
loadings(pca_dm)
## 
## Loadings:
##                         Comp.1 Comp.2 Comp.3
## avg_hindfoot_summer      0.656         0.748
## avg_skull_length_summer  0.581  0.567 -0.585
## avg_skull_width_summer   0.482 -0.818 -0.314
## 
##                Comp.1 Comp.2 Comp.3
## SS loadings     1.000  1.000  1.000
## Proportion Var  0.333  0.333  0.333
## Cumulative Var  0.333  0.667  1.000
summary(pca_dm)
## Importance of components:
##                           Comp.1    Comp.2    Comp.3
## Standard deviation     1.1650734 0.9573210 0.8521387
## Proportion of Variance 0.4524654 0.3054879 0.2420468
## Cumulative Proportion  0.4524654 0.7579532 1.0000000
screeplot(pca_dm)

eigenvals(pca_dm) 
##    Comp.1    Comp.2    Comp.3 
## 1.3573961 0.9164636 0.7261403
body_size_dm <- (pca_dm$scores[,1])

Data exploration

plot(avg_weight_sampling_period ~ body_size_dm, 
     data =dm)

Extracting residuals

bodycond_correlation_dm <- lm(dm$avg_weight_sampling_period ~ body_size_dm)
dm$bodycondition <- bodycond_correlation_dm$residuals
hist(dm$bodycondition)

priors

get_DMbodycondition_priors  <- get_prior(bodycondition ~ infection+sex*scale(julian_date)+ (1|tag_left), 
                      data=dm,
                      family=gaussian()); get_DMbodycondition_priors 
DMbodycondition_priors <- c(set_prior("normal(0,5)", class ="Intercept"),
                             set_prior("normal(0,5)", class="b")) 

model formula

DMbodycondition_formula <- bf(bodycondition ~ infection+sex*
                                 scale(julian_date, center=TRUE, scale=TRUE)+(1|tag_left), 
                               family =gaussian())

Fit model (work in progress)

DMbodycondition <- brm(DMbodycondition_formula, 
                       data =dm, 
                       family=bernoulli(), 
                       prior =DMbodycondition_priors,
                       iter =2000,
                       warmup =1000, 
                       seed=123,
                       thin=2,
                       chains = 4,
                       cores = 2, 
                       save_pars = save_pars(all = TRUE), 
                       sample_prior = "yes")
saveRDS(DMbodycondition, "models/DMbodycondition_2025.RDS")

model validations (convergence)

stan_trace(DMbodycondition$fit)
## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(DMbodycondition$fit) 
## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(DMbodycondition$fit) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(DMbodycondition$fit)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(DMbodycondition$fit, separate_chains = TRUE) 
## 'pars' not specified. Showing first 10 parameters by default.

model investigation

summary
summary(DMbodycondition)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: bodycondition ~ infection + sex * scale(julian_date, center = TRUE, scale = TRUE) + (1 | tag_left) 
##    Data: dm (Number of observations: 184) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 2;
##          total post-warmup draws = 2000
## 
## Group-Level Effects: 
## ~tag_left (Number of levels: 88) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.97      0.18     0.63     1.32 1.00     1063     1608
## 
## Population-Level Effects: 
##                                              Estimate Est.Error l-95% CI
## Intercept                                        0.06      0.35    -0.63
## infectionParasitized_Before                     -0.17      0.34    -0.83
## infectionUninfected_Before                       0.07      0.35    -0.63
## sexM                                            -0.02      0.32    -0.66
## scalejulian_datecenterEQTRUEscaleEQTRUE         -0.71      0.21    -1.11
## sexM:scalejulian_datecenterEQTRUEscaleEQTRUE     0.11      0.25    -0.37
##                                              u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept                                        0.74 1.00     1608     1644
## infectionParasitized_Before                      0.48 1.00     1747     1929
## infectionUninfected_Before                       0.76 1.00     1519     1726
## sexM                                             0.60 1.00     1861     1658
## scalejulian_datecenterEQTRUEscaleEQTRUE         -0.30 1.00     1877     1788
## sexM:scalejulian_datecenterEQTRUEscaleEQTRUE     0.59 1.00     1817     1771
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.30      0.09     1.13     1.50 1.00     1354     1607
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
R-squared
DMbodycondition |> bayes_R2(summary =FALSE) |> median_hdci()

post-hoc analysis

Infection
posteriors <- as.data.frame(DMbodycondition) 
covariate_samples <- posteriors |> select(b_Intercept, b_infectionParasitized_Before,b_infectionUninfected_Before) 
combined_effect <- rowSums(covariate_samples) 
hdi(combined_effect, credMass =0.95)
##      lower      upper 
## -0.9724559  0.7825951 
## attr(,"credMass")
## [1] 0.95
DMbodycondition |> emmeans(~infection, type = "response") 
##  infection           emmean lower.HPD upper.HPD
##  Parasitized_After   0.0441    -0.513     0.671
##  Parasitized_Before -0.1066    -0.737     0.444
##  Uninfected_Before   0.1209    -0.263     0.506
## 
## Results are averaged over the levels of: sex 
## Point estimate displayed: median 
## HPD interval probability: 0.95
DMbodycondition |> emmeans(~infection, type = "response") |> 
  pairs() |> summary(by = NULL, adjust = "tukey", infer=TRUE) 
sex
DMbodycondition |> emmeans(~sex, type = "response")
## NOTE: Results may be misleading due to involvement in interactions
##  sex  emmean lower.HPD upper.HPD
##  F   0.03317    -0.467     0.545
##  M   0.00735    -0.407     0.426
## 
## Results are averaged over the levels of: infection 
## Point estimate displayed: median 
## HPD interval probability: 0.95
DMbodycondition |> emmeans(~sex, type = "response") |> 
  pairs() |>
  summary(by = NULL, 
          adjust = "tukey", 
          infer=TRUE)  
## NOTE: Results may be misleading due to involvement in interactions
julian date
DMbodycondition |> emtrends(var = "julian_date", type = "response") |> pairs(by='infection') |>
  summary(by = NULL, adjust = "tukey", infer=TRUE) 

partial effects plots

DMbodycondition |> conditional_effects(spaghetti = TRUE, ndraws =200)

summary figure

Red backed voles
rbv.infect <- as.data.frame(emmeans(RBVbodycondition, ~infection)) 

rbv.infect.obs <- drop_na(rbv, infection) |> 
  mutate(Pred =predict(RBVbodycondition, re.form =NULL, type ='response'), 
         Resid =residuals(RBVbodycondition, type ="ordinary"), 
         Fit =Pred+Resid) 
rbv.bodycondition.plot <- ggplot(data =rbv.infect, aes(x=infection, y=emmean)) + 
  geom_jitter(data =rbv.infect.obs, aes(x= infection, y =Fit[,1]), 
              width=0.02,
              color ="#7E4640", 
              alpha =0.15) +
  geom_pointrange(aes(ymin =lower.HPD, 
                      ymax =upper.HPD), 
                      color ="black")+ 
  xlab("") + ylab("Body condition") + 
  scale_y_continuous(limits =c(-10,10)) +
  scale_x_discrete(limits =c("Uninfected_Before","Parasitized_Before","Parasitized_After")) +
  theme_classic() +
  theme(legend.position = "none", 
        axis.text.x = element_blank()); rbv.bodycondition.plot

Deer mice
dm.infect <- as.data.frame(emmeans(DMbodycondition, ~infection)) 

dm.infect.obs <- drop_na(dm, infection) |> 
  mutate(Pred =predict(DMbodycondition, re.form =NULL, type ='response'), 
         Resid =residuals(DMbodycondition, type ="ordinary"), 
         Fit =Pred+Resid) 
dm.bodycondition.plot <- ggplot(data =rbv.infect, aes(x=infection, y=emmean)) + 
  geom_jitter(data =rbv.infect.obs, aes(x= infection, y =Fit[,1]), 
              width=0.02,
              color ="#7C8081", 
              alpha =0.25) +
  geom_pointrange(aes(ymin =lower.HPD, 
                      ymax =upper.HPD), 
                      color ="black") + 
  xlab("Infection status") + ylab("Body condition") + 
  scale_y_continuous(limits =c(-10,10)) +
  scale_x_discrete(limits =c("Uninfected_Before","Parasitized_Before","Parasitized_After")) +
  theme_classic() +
  theme(legend.position = "none"); dm.bodycondition.plot

merge plots
rbv.dm.bodycondition.plot <- (rbv.bodycondition.plot / dm.bodycondition.plot) + plot_annotation(tag_levels = 'A') ; rbv.dm.bodycondition.plot

save plot
ggsave("outputs/2025_bodycondition2.pdf", width =5, height =5, units ="in", dpi =300) 
rbv.dm.bodycondition.plot 
dev.off()